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Abstract 


The  advantages  of  solving  two-dimensional  potential  problems  using  spectral 
boundwy  integral  methods  are  examined.  Using  fast  Fourier  transforms,  we  ex¬ 
pand  the  spatizd  coordinates  x  and  y  using  an  arclength  par2uneter  s.  This  spectral 
representation  is  very  accurate  when  the  geometry  is  smooth  and  nodal  spacing 
is  uniform.  Two  spectral  formulations  are  outlined.  One  is  based  on  Baker’s  [11] 
integration  scheme  at  every  other  node  to  avoid  the  kernel  singularities,  and  the 
other  is  based  on  the  kernel  desingularization  of  Roberts  (12).  An  error  analysis 
and  convergence  studies  for  several  j^eometries  are  shown. 
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1  Introduction 


Boundary  integral  methods  and  spectral  methods  are  two  of  the  most  powerful 
techniques  in  computational  mechanics.  In  this  paper  we  combine  these  approaches 
to  obtain  a  very  accurate  algorithm  for  potential  problems.  The  governing  equation 
for  potential  problems  is  elliptic,  requiring  one  condition  on  the  enclosing  bound¬ 
ary.  The  boundary  integral  method  only  solves  for  an  unknown  on  the  boundary, 
which  reduces  the  effective  dimension  of  the  potential  problem  by  one  and  signif¬ 
icantly  decreases  the  computational  effort  unless  the  contour  is  highly  distorted. 
In  addition,  if  we  use  spectral  representation  for  the  boundary  values,  we  may  get 
exponentially  accurate  solutions. 

The  solutions  of  two-dimensional  potential  problems  can  be  described  by  a  com¬ 
plex  analytic  function.  The  boundary  integral  method  can  then  be  derived  from 
Cauchy’s  integral  theorem  (Greenhow,  Vinje,  Brevig,  and  Taylor  [1],  Vinje  and 
Brevig  [2],  Lai  2uid  Hromadka  [3]),  rather  than  the  more  general  method  using  a 
free-space  Green’s  function  (Longuet-Higgins  and  Cokelet  [4]).  Applying  the  real 
and  imaginary  parts  of  the  Cauchy  integral  theorem  results  in  two  real  algebraic 
systems.  While  Cauchy’s  integral  and  the  Green’s  function  formulations  are  not 
directly  compared  here.  Bold  and  Peregrine  [5,6]  indicate  that  the  former  method 
[1,2,3]  is  clearly  superior.  Schultz  [7,8]  shows  that  the  error  of  the  complex  formu¬ 
lation  can  be  reduced  by  solving  these  two  systems  of  equations  in  a  least-square 
sense  for  nonlinear  breaking  wave  simulations.  Schultz  and  Hong  [9]  compare  the 
results  for  the  overdetermined,  strong,  and  weak  systems  using  a  piecewise-linear 
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method.  They  concluded  that  the  overdetermined  system  is  superior  to  others  in 
many  cases. 

In  general,  an  integral  equation  is  solved  by  a  numerical  model  that  assumes 
the  boundary  is  composed  of  piecewise-polynomial  curves  (panels)  and  the  known 
and  unknown  boundary  values  are  approximated  as  piecewise-continuous  functions 
along  the  boundary.  Most  complex  variable  boundary  integral  methods  adopt 
piecewise-linear  representations  of  the  a>mplex  functions.  A  linear  interpolation 
then  results  in  second-order  accurate  integration  and  typically  gives  second-order 
accurate  solutions  for  the  boundary  integral  solution  as  well.  To  improve  accuracy, 
we  expand  the  known  and  unknown  boundary  vriues  and  the  boundary  shape  in  a 
global  Fourier  series  using  an  arclength  parameter.  Spectral  methods  then  result 
in  exponentially  accurate  integration  and  give  exponentially  accurate  solutions. 

In  this  paper,  two  spectral  formulations  are  outlined.  One  is  based  on  Baker’s 
[11]  integration  at  every  other  node  point  to  avoid  the  kernel  singularities,  and  the 
other  is  based  on  the  kernel  desingularization  of  Roberts  [12].  The  advantages  of 
the  spectral  boundary  int^al  method  are  discussed,  and  the  error  is  analyzed  and 
compared  to  piecewise-linear  method  for  several  test  cases.  The  test  cases  show 
the  effects  of  geometric  curvature,  nodal  spacing,  and  local  solution  gradients  on 
the  solution  error. 
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2  Problem  Formulation 


Figure  1  shows  the  problem  domain  with  the  bounding  surface  dR.  Solutions 
of  potential  problems  ^  and  tj;  can  be  written  as 

0{2)  =  <f>+  .  (1) 

Here,  z  =  x  +  iy  represents  the  two  spatial  coordinates.  Cauchy’s  theorem  then 
gives 

where  a  is  0  or  2ir  if  the  location  of  the  kernel  singularity,  (k,  is  outside  or  inside 
the  boundary,  respectively.  If  the  kernel  singularity  is  on  the  boundairy  (C*  6  dR), 
a  is  the  included  angle,  and  (2)  is  treated  as  a  principal  value  integral. 

Our  goal  is  to  convert  (2)  into  an  algebraic  system  as  in  Schultz  and  Hong  [9]. 
The  algebraic  system  is  formed  by  discretizing  the  integral  and  letting  the  kernel 
singularity  approach  each  of  the  N  nodal  points,  —*■  ^k-  A  speciad  limiting  pro¬ 
cess  is  then  needed  to  evaluate  the  integration  near  Qk- 

2.1  Spectral  Method  I 

Following  the  method  outlined  by  Baker  [11],  exponential  convergence  can  be 
achieved  by  evaluating  the  principal  value  int^rad  of  (2)  using  the  trapezoidal 
rule.  To  obtain  this  convergence,  the  integrad  must  be  periodic  and  integrated 
with  respect  to  a  real  variable,  in  this  case,  an  au-clength  pairameter,  s.  The  integral 


becomes 


Baker  [11]  indicates  this  expression  will  give  exponentially  convergent  results  when 
trapezoidal  integration  is  used  for  a  smooth,  periodic  integraind  and  the  shape  (and 
hence  dz/ds)  is  known  exactly.  In  contour  d)rnamics  problems,  for  example,  the 
bounding  surface  dR  must  be  calculated  from  evolution  equations.  However,  if 
the  contour  is  known  to  be  smooth,  z(s)  and  hence  dzfds  can  be  determined  with 
exponential  accuracy  if  spectral  methods  are  used.  We  have  performed  these  cal¬ 
culations  using  a  standard  fast  Fourier  trsmsform  (FFT)  routine  for  both  x(s)  and 
y(s)  to  obtain  exponential  accuracy. 


In  (3),  the  integrand  is  singular  as  z  approaches  The  singular  part  of  the  in¬ 
tegral  is  avoided  to  smooth  the  integrand  by  evaluating  the  principal- value  integral 
at  every  other  point.  The  algebraic  equations  can  be  written  as 


where 


=  0 

i=i 


for  A  =  1 . . . ,  TV  , 


if|i-*|isodd 


r,*4  —  <  —in 


if  j  =  k 


[  0  if  \j  —  l:j  is  even  and  j  ^  k  . 

If  ^  is  given  on  the  boundary,  tp  can  be  calculated  from  (4)  and  (5).  Since  the 
integration  is  carried  out  at  every  other  point,  unlike  the  piecewise-linear  repre¬ 
sentation  of  [9],  many  Tji,  are  zero  (nearly  one-half).  Obviously,  this  saves  much 
time  (nearly  one-half)  in  solving  the  algebraic  system. 
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2.2  Spectral  Method  II 


Roberts  [12]  used  a  desingularized  kernel  in  his  vortex  formulation.  Generally,  it 
is  difficult  to  find  a  suitable  desingularized  form  of  a  kernel  in  an  integral  equation, 
but  in  the  complex  formulation  it  is  relatively  simple.  The  Cauchy  integral  equation 
(3)  can  be  easily  rewritten  as 

J  z(s)  —  Zk  as 

where  the  principal  value  integral  can  be  replaced  by  the  closed  contour  integral 
since  the  integrand  is  no  longer  singular.  When  z  approaches  z^,  the  integrand 
approaches  d^Jds  at  the  ibth  node.  Therefore,  this  kernel  does  not  exhibit  singular 
behavior  as  z  approaches  z*.  Integral  equation  (6)  is  converted  to  the  following 
sets  of  equations  for  A:  =  1  ...N: 


N 


^j-Pkfdz, 


T-rT:(x)i  +  =  0 


for  A:  =  1, . . . ,  JV  , 


(7) 


where  N  is  the  number  of  nodes.  The  desingularized  algebraic  system  (7)  requires 
the  derivative  of  and  hence  effectively  becomes  a  differential  system.  To  evaluate 
these  derivatives  spectrally,  we  use  a  trigonometric  interpolation  of  ^  [13]: 


N 


where 


i=i 


Ci  =  ^  sin  7r(s  -  Sj)  cot  ^(s  -  s^)  , 


(8) 


(9) 


and  the  derivative  of  Cj  is 


ds 


cot  ^{s<  -  Sj)  if  i  j 

0  if  *  =  3 


(10) 
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Then,  (7)  becomes 


N 


EfikA  =  0 


for  A  =  , 


(11) 


where  the  influence  coeflScients  Fjib  are  now 

I  _ ar  •  *  •  . 


if  j  =  k  . 


(12) 


Unlike  spectral  method  I,  the  desingularized  kernel  is  evaluated  at  every  nodal 
point,  and  hence  the  matrix  is  full. 


3  Numerical  Investigation 

Comparisons  of  accuracy,  convergence  rate,  and  computing  time  are  made  for 
spectral  method  I,  spectral  method  II,  and  the  piecewise-lineax  method.  The  ef¬ 
fects  of  nonimiform  nodal  spacing,  contour  shape,  and  nearly  singular  solution 
characteristics  on  the  solution  error  are  investigated. 

To  solve  the  matrix  problem,  a  conjugate  gradient  iterative  technique  is  used. 
This  technique  solves  the  problem  in  an  order  of  UVL  operations,  where  U  is  the 
ntimber  of  unknowns,  V  is  the  number  of  equations,  and  L  is  the  number  of  itera¬ 
tions.  The  computational  savings  of  the  iterative  technique  would  be  important  if 
time  marching  were  desired,  especially  since  a  good  initial  guess  is  available  from 
the  previous  time  step.  However,  even  with  a  homogeneous  initial  guess,  a  typical 
solution  for  a  nonsingular  matrix  requires  less  than  10  iterations. 
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According  to  Schultz  and  Hong  [9],  the  overdetermined  system  for  the  piecewise- 
linear  method  is  generally  most  accurate,  and  the  convergence  rate  of  the  overdeter¬ 
mined  system  can  be  twice  that  of  the  weak  or  strong  algebraic  systems.  Similarly, 
the  real  and  imaginary  parts  of  (5)  for  spectral  I  and  (12)  for  spectraJ  II  can  be 
imposed  as  well.  Although  we  find  that  for  uniformly  spaced  nodes  the  conver¬ 
gence  rates  of  both  spectral  methods  are  nearly  identical  for  all  algebraic  systems, 
the  overdetermined  system  gives  more  stable  results  for  distorted  contours  and 
nonuniform  nodal  spacing  [14].  We  use  an  overdetermined  system  for  all  methods 
unless  otherwise  indicated. 

The  test  function  is  ^  =  sin  z  unless  otherwise  stated.  The  real  part  of  i.e., 
4>,  is  prescribed  at  each  node;  ip  is  then  determined  from  the  boundary  integral 
method  and  compared  to  the  known  We  perform  all  calculations  on  an 

Apollo  DN4000  Workstation  using  double  precision  (16  digits). 

3.1  Rate  of  Convergence  and  Computing  Time 

The  root  mean  square  errors,  £^2,  evaluated  at  the  nodal  points  and  the  cor¬ 
responding  computing  times  are  examined  as  a  function  of  the  number  of  nodes 
to  determine  convergence  rates  and  computational  efficiency.  Figure  2  compares 
£2  for  imiformly  spaced  nodes  on  a  circular  contour  for  three  methods:  spectral 
I,  spectral  II,  and  piecewise-Iinear.  The  spectral  methods  are  much  more  accurate 
than  the  piecewise-Iinear  method,  and  spectral  II  is  slightly  better  than  spectral  I. 


The  effect  of  the  contour  shape  on  the  efficiency  is  examined  for  elliptic  con¬ 
tours  with  vairious  eccentricities  (c  =  0.0  to  0.995).  The  diameter  of  the  major 
zacis  is  fixed  at  2.0  and  centered  on  the  x  axis.  Figure  3  shows  solutions  for  var¬ 
ious  elliptic  contours  when  N  is  16.  Here,  the  arclength  parameter  spacing  is 
determined  by  As  =  rAO,  where  A6  is  constant.  As  the  eccentricity  increases, 
the  Ei  of  spectral  I  is  nearly  constant  up  to  e  =  0.7  and  then  increases  rapidly; 
that  of  spectral  11  also  remains  nearly  constant  for  small  e  but  decreases  at  e  =  0.7. 

The  nodal  spacing  s  varies  as  the  eccentricity  increases,  and  the  integration 
neair  the  singular  point  is  not  accurate  using  spectral  I,  which  explains  the  increase 
in  Ej  after  e  =  0.7.  On  the  other  hand,  the  integrand  in  spectral  I  has  a  limiting 
form  as  z  goes  to  z*.  Therefore,  Ej  is  independent  of  eccentricity.  Both  spectral 
methods  are  much  more  accurate  than  the  piecewise-linear  method  except  near  e 
=  1 ,  where  the  piecewise-linear  method  is  more  accurate  than  spectral  I. 

In  Figure  4,  error  vs.  CPU  time  is  shown  for  an  elliptic  contour  (c  =  0.866). 
Generally  the  CPU  time  of  spectral  I  is  the  smallest  for  a  given  N,  although  the 
greater  accuracy  of  spectral  II  makes  it  more  computationaJly  efficient.  Figure  4 
clearly  shows  that  both  spectral  methods  are  more  efficient  than  the  piecewise- 
linear  method. 
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3.2  Nonuniform  Nodal  Spacing 


Sometimes  nonuniform  nodal  spacing  is  desirable,  as  when  solution  refinement 
is  needed  due  to  large  solution  gradients.  The  nodal  spacing  also  varies  in  problems 
with  convective  nodes  such  as  free  surface  flow  problems.  This  can  cause  numerical 
instabilities,  requiring  a  filtering  or  a  regriding  scheme  to  maintain  uniform  spacing. 

In  the  previous  section,  we  note  that  as  e  increases,  the  nodal  spaunng  about 
an  arclength  varies,  but  the  accuracy  remains  good  for  spectral  II.  In  those  cases, 
even  though  the  nodal  spacing  about  the  arclength  is  nonuniform,  its  variation 
is  very  smooth  and  the  nodal  spacing  about  6  is  uniform.  However,  if  only  one 
node  is  moved,  the  abrupt  change  in  nodal  spacing  seriously  affects  the  fast  Fourier 
transform  accuracy.  Here,  two  cases  are  considered:  one  when  an  abrupt  chzmge 
occurs  in  nodal  spacing,  and  the  other  when  the  nodal  spacing  varies  slowly. 

The  first.case  is  shown  in  Figure  5  for  a  circular  contour.  When  node  1  is  moved 
a  fraction  p  of  the  uniform  nodal  spacing  A9  =  2v/N  towards  node  2,  then  the 
accuracy  decreases  rapidly  as  p  increases  for  all  methods.  Figure  5  compares  the  re¬ 
sults  of  the  three  methods  when  p  is  0.1.  For  N  >  16,  the  piecewise-linear  method 
gives  the  most  acairate  solution,  and  for  all  N,  spectral  II  gives  more  accurate 
solutions  than  spectral  I,  as  expected.  Figure  6  compares  the  results  for  various  p 
when  N  is  16.  Near  p  =  0,  the  spectral  methods  give  more  accurate  results  than 
the  piecewise-linear  method,  and  the  accuracy  increases  rapidly  for  0  <  p  <  0.1. 
But  as  p  increases,  the  spectral  methods  become  inaccurate.  The  accuracy  of  the 
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piecewise- linear  method  does  not  degrade  as  significantly  as  p  increases. 


We  conclude  that  the  piecewise-linear  method  is  less  succeptible  to  an  abrupt 
change  of  nodal  spacing.  The  inaccuracies  of  the  spectral  methods  are  due  to 
integration  errors  as  well  as  FFT  routine  errors.  Integration  errors  affect  the  accu¬ 
racy  of  spectral  I  more  significantly  than  spectral  II  because  singular  contribution 
on  the  integral  doesn’t  cancel  for  uneven  nodal  spacing  in  spectral  I,  while  FFT 
routine  errors  affect  the  accuracy  of  both  spectral  methods.  Hence  spectral  II  can 
tolerate  larger  nodal  spacing  irregularities  than  spectral  I. 

The  second  case  is  shown  in  Figure  7,  which  compares  the  results  of  the  three 
methods  for  an  elliptic  contour  of  slowly  varying  nodal  spacing.  Here,  AB  varies 
as  a  cosine  function,  b  is  the  amplitude  of  the  cosine  function,  and  N  =  32.  Both 
spectral  methods  are  more  accurate  than  the  piecewise-linear  method,  and  spectral 
method  II  is  significantly  better  than  spectral  I. 


3.3  Distorted  Contours 

The  computational  domain  can  become  very  distorted.  For  example,  in  wave 
modulation  problems,  long  waves  are  combined  with  short  waves.  After  conformal 
mapping,  a  contour  similar  to  the  following  is  obtsuned: 

X  =  r(^)  cos  9  , 

y  =  r(^)  sin  9  , 
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r(0)  —  1  +  a  cos  k6  . 


In  Figure  8,  the  domain  geometry  is  shown  when  k  is  8.  Figure  9  compares  the  re¬ 
sults  of  the  three  methods.  In  this  case,  spectrcil  II  gives  the  most  accurate  results, 
and  when  N  >  128,  spectral  I  is  more  accurate  thzm  the  piecewise-linear  method. 
Figures  10  and  11  show  the  results  for  different  k  using  spectral  I  and  spectral 
II,  respectively.  As  expected,  as  k  increases,  the  accuracy  decreases,  because  the 
domain  geometry  is  more  contorted.  Figure  12  compares  the  three  methods  for 
vMying  a  when  A  is  8  and  N  is  64.  When  a  =  0  (circular  contour),  both  spectral 
methods  are  more  accurate  than  the  piecewise-linear  method,  as  shown  previously. 
As  a  increases,  in  spectral  I  increases  rapidly,  and  when  a  >  0.1,  the  accuracy 
of  spectral  I  is  worse  than  that  of  the  piecewise-linear  method.  In  spectral  II,  the 
Fj  also  increases,  but  the  accuracy  is  still  better  than  the  piecewise-linear  method. 
From  these  test  cases,  spectral  II  is  shown  to  be  the  most  efficient  method  for 
distorted  contours. 

The  boundary  integral  accuracy  is  known  to  deteriorate  in  the  vicinity  of  con- 
toTir  comers  due  to  the  infinite  geometric  curvature  (Schultz  and  Hong  [9]).  With 
spectral  methods,  since  the  FFT  for  x(s)  and  y(s)  converges  very  slowly,  our  re¬ 
sults  show  very  poor  results.  Therefore,  we  conclude  that  spectral  methods  offer 
no  advantages  for  contours  with  comers.  The  results  are  not  shown  here. 
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3.4  Desingularization 


We  use  the  same  strategy  for  the  locating  the  kernel  singularities  outside  the 
contour  as  Schultz  and  Hong  [9].  The  strategy  uses  the  perpendicular  bisector 
of  the  straight  contour  between  adjacent  nodes.  This  method  has  two  potential 
difficulties,  as  pointed  out  by  Schultz  and  Hong  [9]:  if  /  is  too  small,  the  singular 
point  may  lie  inside  a  convex  contour,  but  if  /  is  too  large,  the  singular  point  may 
lie  inside  the  domain  of  a  highly  contorted  contour. 

Figure  13  shows  the  error  computed  for  several  values  of  /  in  spectral  I,  where 
the  nodes  are  placed  evenly  on  a  circrilar  contour  with  ^  =  sin  z.  Figure  14  shows 
the  error  computed  for  various  values  of  /  in  spectral  I  where  N  is  16.  In  this  case, 
since  the  desingularized  kernel  is  used,  integration  is  carried  out  at  every  point. 
Generally,  the  algebraic  system  becomes  less  diagonally  dominant  as  /  increases, 
requiring  more  iterations  and  a  tighter  tolerance  on  the  convergence  parameter  to 
zxhieve  the  desired  accuracy. 


4  Conclusions 

Our  numerical  investigation  shows  that  the  spectr2tl  II  is  the  most  efficient  amd 
accurate  in  almost  all  cases.  For  smooth  contours  with  uniform  nodal  sparing, 
both  spectral  methods  are  more  accurate  thain  the  piecewise-linear  method.  The 
piecewise-linear  method  is  less  succeptible  to  an  abrupt  change  in  the  nodal  spac- 
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ing  than  spectral  methods,  but  when  the  nodal  spacing  varies  slowly,  spectral  II 
is  most  accurate.  For  distorted  contours,  spectral  II  is  the  most  efficient  method, 
because  it  gives  more  accurate  results  for  the  same  N.  Both  spectral  methods 
have  difficulties  handling  contours  with  comers.  As  for  efficiency,  both  spectral 
methods  need  less  computing  time  than  the  piecewise-linear  method,  and  spectral 
II  is  the  fastest  scheme.  The  accuracy  of  spectral  II  may  be  further  improved  by 
filtering  or  a  node  redistribution  scheme,  especially  when  discontinuities  in  nodal 
spacing  or  contour  shape  are  encountered. 

The  spectral  approach  may  be  extended  to  periodic  three-dimensional  potentiad 
problems  using  integral  equations  derived  from  Green’s  theorem  or  the  normal 
derivative  of  Green’s  theorem.  Hence,  the  spectral  boundary  integral  method  may 
yield  greater  efficiency  and  accuracy  for  three-dimensional  problems  as  well. 
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LIST  OF  SYMBOLS 


a  :  amplitude  of  cosine  variation  of  r(d)  in  complicated  contour 

b  :  amplitude  of  cosine  variation  of  Ad 

:  root  mean  square  error 

c  :  eccentricity  of  an  ellipse,  y/A^  —  B^/A^  A  =  long  axis,  B  =  short  axis 

/  :  desingularization  parameter  for  locating  kernel  singular  point 

Ijit  :  discretized  kernel 

k  :  niunber  of  bumps  in  complicated  contour 

N  :  total  number  of  nodes 

p  :  amounts  of  uniform  AO  moving  from  node  1  to  node  2 

R  :  domain  of  the  problem 

3  :  curvUinear  coordinate  on  the  contour 

U,V  :  row  and  column  number  of  an  overdetermined  matrix 

z(s),y(s)  :  X  and  y  coordinate  of  the  contour 

z  :  complex  coordinate  of  a  point 

Zj  :  j**  node 

Greek  Letters 

a  :  a  constant  (alpha) 

:  complex  potential  (beta) 

Tjk  :  integral  equation  influence  coefficient  (upper  case  Gamma) 

Ck  :  control  point  (zeta) 
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:  tangent  angle  of  the  contour  with  positive  x-axis  (theta) 
<l>  :  real  part  of  0  (phi) 

Ip  :  imaginary  part  of  0  (psi) 

dR  :  contour  of  computing  domain 
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